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ABSTRACT 

We present a new method for analyzing the global stability of the Sedov-von Neumann- Taylor self- 
similar solutions, describing the asymptotic behavior of spherical decelerating shock waves, expanding 
into ideal gas with density oc r~". Our method allows to overcome the difficulties associated with 
the non-physical divergences of the solutions at the origin. We show that while the growth rates of 
global modes derived by previous analyses are accurate in the large wave number (small wavelength) 
limit, they do not correctly describe the small wave number behavior for small values of the adiabatic 
index 7. Our method furthermore allows to analyze the stability properties of the flow at early times, 
when the flow deviates significantly from the asymptotic self-similar behavior. We find that at this 
stage the perturbation growth rates are larger than those obtained for unstable asymptotic solutions 
at similar [7,w]. Our results reduce the discrepancy that exists between theoretical predictions and 
experimental results. 

Subject headings: hydrodynamics — instabilities — shock waves 



1. INTRODUCTION 

Shock waves play a crucial role in the evolution of a 
wide variety of astrophysical systems, such as the inter- 
stellar medium and the inter- galact ic medium (for review 
see, e.g.. lOstriker fc McKee Ill988j) . The study of shock 
wave properties, and in particular of shock wave stabil- 
ity, is therefore relevant for a wide variety of astrophys- 
ical phenomena. The stability of planar steady shocks 
pr opagating in t o a u niform medium was demonstrated 
bv'Eripenbeckl l|1962D . The stability of ID shocks (with 
planar/cylindrical/spherical symmetry) propagating into 
both uniform and non-uniform media is commonly stud- 
ied by analyzing self-similar shock solutions of the hy- 
drodynamic equations, since self-similar solutions often 
allow analytic stability analysis and since in many cases 
self-similar solution s describe the asymptotic behavior of 
the flow (see, e.g., Zcldovic h fc R,aizer 1119681) . 

Ryu & Vishniac (1987,1991) have studied the stability 
of the self-similar Sedov-von Neumann- Taylor solu tions 
llSedovllT94^ iVon Neum^^ 179471: iTavlor llT9Rfil) . de- 
scribing spherical shocks propagating into an ideal gas 
with a power-law density profile p oc r~". They have 
studied both the "local" and the "global" stability of 
the solutions. The term "local" stability analysis refers 
to approximate analysis of stability where flow prop- 
erties are considered in a small region, where the spa- 
tial dependence of flow variables may be approximated 
as linear, and boundary conditions are neglected. The 
term "global" stability refers to stability against pertur- 
bations with time independent spatial structure, which 
satisfy the flow boundary conditions. Analysis of sta- 
bility against global modes of perturbation is often not 
useful for general time dependent flows, with time de- 
pendent spatial structure. Self-similar flows are an ex- 
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ception, since the spatial structure of such flows is time 
independent. 

The amplitude of global modes of perturbations to the 
Sedov-von Neumann- Taylor solutions evolves with time 
as i*. For expansion into uniform media (w = 0), the 
Sedov-von Neumann- Taylor shocks were found to be sta- 
ble (i?e(s) < 0) for gas adiabatic index 7 > 1.2, and 
unstable (for a limited range of perturbation wave num- 
bers) for 7 < 1.2, where the high density region behind 
the shock is sufficiently thin. Qualitatively similar results 
are obtained from the analysis of Ryu & Vishniac for uj 
values in the range < w < (7 — 7)7(7 -I- 1). For larger 
values of lo, (7 — 7)/(7 + 1) < w < 3, the solutions are 
"hollow," i.e., the density vanishes at some finite distance 
(behind the shock) away from the origin. The stability 
of "hollow" solutions was examined by Goodman (1990), 
who found that the flow is unstable for large perturba- 
tion wave numbers Z, with gro wth rate scaling a s V-^"^ . 
The local stability analysis of 'Rvu fc Vishniac I l)1991f) 
showed that the flow in the non-hollow solutions is con- 
vectively unstable (to local perturbations) in the range 
3/7 < a; < (7 — 7)7(7-1-1). The growth rate of the per- 
turbation was studied nu merically in one speci f ic cas e 
(a; = 0,7 = 1.1,1 = 40) bv lMac Low fc Norman I l|l99^. 
and was found to be in good agreement with the predic- 
tion of Ryu & Vishniac. 

Shock waves propagating into a steep density gra- 
dient > 3, are accelerating. The asymptotic 
behavior of such shocks is not described by the 
Sedov-von Neumann- Taylor solutions, but rather by 
a different family o f self-similar solutions, derived by 
iWaxman k. Shvartsi l(T993i'l . The stability of these so- 
lutions ag ainst global modes of perturbati ons was ex- 
amined by 'Sar i. Waxman fc Shvartsi l|200 n1. who found 
that shock waves that accelerate at a rate faster than the 
critical rate RR/R^ — 1 are unstable for small and inter- 
mediate wave numbers, whereas shock waves accelerating 
at a slower rate, RR/R^ < 1, are stable for most wave 
numbers (here, R{t) denotes the time dependent shock 
radius). They have further shown that perturbations of 
small wave numbers grow or decay monotonically, while 
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perturbations of intermediate and high wave numbers os- 
cillate in time. 

The stability of spherical shocks propagating into 
a uniform mediu m was studied experimentally by 
iGrun et al. 1 ^91). For shocks propagating in Xenon, 
with adiabatic index 7 = 1.06±0.02, perturbations grow- 
ing as a power low of time were found, while in Nitrogen, 
7 = 1.3 ± 0.1, shocks were found to be stable. These 
results are in qualitativel y agree ment with the predic- 
tions of ' Rvu fc Vishniac I (|1987). However, the maxi- 
mum growth rate measured was higher then the Ryu & 
Vishniac prediction by a factor of ^ 1.7, and was ob- 
tained at lower than predicted wave numbers, I 10 
instead of Z ~ 60. 

The cause for the discrepancy between the experimen- 
tal results and the theoretical predictions is yet unknown. 
On the one hand, there are some preliminary claims 
that the the expe rimental results of Gr un et al. (1991) 
are not accurate ijHansen et al. II2OO2D . On the other 
hand, there are some difficulties in the theoretical anal- 
ysis. Briefly, the Sedov-von Neumann- Taylor solutions 
for cj < (7 — 7)7(7 + 1) show non-physical divergences 
near the origin, reflecting the fact that for any physi- 
cal initial conditions the flow deviates significantly from 
the asymptotic self-similar solution at all times for suffi- 
ciently small radii, and the analysis of Ryu & Vishniac 
(1987,1991) relies on applying boundary conditions for 
the perturbation equations at the origin. This leads, for 
example, to non physical results of the stability analysis 
for small wave number (large wave length) perturbations 
(see §|21l. 

In this paper we present a new method, free of such dif- 
ficulties, for determining shock stability against global 
perturbations. We analyze the stability of solutions 
which coincide with the Sedov-Taylor-von Neumann so- 
lutions everywhere except for some region of finite mass, 
m, near the origin. We show that the global stability of 
the flow is independent of the details of the solution in the 
region where it deviates from self-similarity, and obtain 
the stability properties of the self-similar solution by tak- 
ing the limit m — > 0. We flnd that while the growth rates 
derived by the method of Ryu & Vishniac approach our 
growth rates for large wave number (small wavelength) 
perturbations, they do not correctly describe the small 
wave number behavior for small values of the adiabatic 
index 7. We furthermore argue that the growth rates ob- 
tained for flnite values of m describe the stability prop- 
erties of the flow at early times, when the flow deviates 
signiflcantly from the asymptotic self-similar behavior. 
We find that at this stage the perturbation growth rates 
are larger than those obtained for unstable asymptotic 
solutions at similar [7,0;]. 

Several attempts to resolve the discrepancy between 
theoretical and experimental results have been made re- 
cently, based on atomic-physics calculations of radiative 
cooling which suggest that the effective adiabatic index 
of the gas relevant for the experiment is l ower than pre- 
viously estimated (Laming fc Grun |I2002L 2003). These 
recent analyses predict larger growth rates, reducing the 
discrepancy with the experimental growth rates. How- 
ever, the most unstable modes are predicted to have 
larger wave numbers, I ~ 200, increasing the discrepancy 
with the experimental results, I 10. In our analysis, 
the large growth rates obtained in the experiment may 



be explained as due to instabilities of the flow at a stage 
where it significantly deviates from self-similarity (finite 
m). The most unstable modes predicted by our linear 
analysis are still characterized by larger wave numbers 
{I ~ 50) than obtained by Grun et al. (1991). However, 
since the growth rates predicted by our analysis for large 
wave numbers, I ~ 50, in the pre-asymptotic stage are 
large, we expect the evolution of perturbations to become 
non-linear at an early stage. The non-linear interaction 
at large wave numbers may lead to perturbation ampli- 
tudes at smaller wave numbers, I ~ 10, which are larger 
than predicted by linear analysis. 

This article is organized as follows. In §|21we provide a 
brief description of the (unperturbed) self-similar shock 
solutions, with particular emphasis on flow properties 
near the origin. The perturbation analysis method is 
described in § 13 The method that allows to overcome 
the difficulties posed by the divergences near the origin 
is described in § 13.21 The stability analysis results are 
described in § ^ A discussion of our results and of their 
implications to the experimental results is given in § [S] 

2. THE UNPERTURBED SHOCK SOLUTIONS 

Consider the strong explosion problem, where a large 
amount of energy is released at the center of a sphere of 
ideal gas with a density profile decreasing with distance 
from the origin as p = Kr~^ . The energy release drives 
a strong outgoing shock wave. We present below a brief 
description of the derivation of the Sedov-von Neumann- 
Taylor solutions, that describe the asymptotic behavior 
of the flow (approached as the shock radius diverges) for 
a; < 3, and point out some of their properties that are 
relevant for our perturbation analysis. 

2.1. The self-similar solutions 

The hydrodynamic equations describing the flow of an 
ideal gas with adiabatic index 7 in a spherically symmet- 
ric geometry are: 

{dt + udr)p + pr^^dr{r^u) = 0, 

p{dt + udr)u + dri-y^^pc^) = 0, 

idt+udr)h-^C^P^-^) = 0, (1) 

where the dependent variables u, c, and p are the fluid 
velocity, sound velocity, and density respectively. The 
pressure is given via equation of state for ideal gas as: 
p = p(? j^. A self-similar solution to the hydrodynamic 
equations (Q) is a solution of the form : 

p{r, t) = BR^GiO, P{r, t) = BR^R^P{0, (2) 
where ^ = r/R{t) is the dimensionless spatial coordinate 
and the length scale R{t), which is c hosen in the present 
case t o be the shock radius, satisfles (^Zeldovich fc Raizerl 
119681 IWaxman fc Shvarts1IT99l : 

^^S^RocR^. (3) 

The quantities G, C, U, and P defined by these equations 
give the spatial dependence of the hydrodynamic quanti- 
ties. The diverging (exploding) solutions of equation Q 
are : 

r A{t~tor, 5<i, 

R{t) - <^ Ae*/-, S - 1, (4) 
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where a = 1/(1 — S). 

The similarity parameter e is determined by the bound- 
ary conditions. These are determined at the shock, £ = 1, 
by th e Rankine-Hugoniot relations ijLandau fc Lifshitz I 
IIQSTT . These relations applied to the self-similar solution 
imply e = —u! and 

[/(!) = 



7 + 1 



G(l): 



V27(7-l) 

7 + 1 '' 
7 + 1 



7 ■ 



1 



(5) 



Substituting equation jSJ into the hydrodynamic equa- 
tions (^3) and using equation the set of par- 
tial differential equations is reduced to a single or- 
dinary differential equation ijZeldovich fc Raizer 1 1196^ 
iMever-ter-Vchn . .1982) 



dU 
dC 



Ai([/,C) 
A2(C/,C)' 



and one quadrature 

dlog£ _ A{U,C) 



dlogC _ A([/,C) 



dU Ai{U,C) dC 
The functions A,Ai, and A2 are 

A = C2 - (1 - C/)2, 



A2(C/,C)- 



(6) 



(7) 



Ai= t/(l - U){1 -U - 



2[(c 



a — 1 , 
a 

l)/a], 



C'^(3C/ 



A2=C[(1-C/)(1-C/-^^) 
^^t/(2(l-C/) + ^)-C^ 



(7- lV + 2[(a- 



a 
I) /a] 



27 \-U" 
and G is given implicitly by 

C-\l - [/)>G^-i+>^3>-2 ^ const. 

with 

A 



7Ct; 



(8) 



(9) 



(10) 



The value of i5 is determined by the conservation of en- 
ergy. Requiring the total energy to be time independent 
gives 5 = (w— 3)/2 (for w > 3 the energy in the Sedov-von 
Neumann- Taylor solutions diverges, and 5 is determined 
by a different argument, see Waxman & Shvarts 199^. 

2.2. The behavior near the origin for u < oJuoUow 

The flow properties are qualitatively different in the 
two regimes w < (7 - 7)/ (7 -I- 1) = Uhoiiow and ujhoUow < 
w < 3. For CO < LUhoiiow U tends to I/7 and C tends 
to infinity as £ tends to zero. For to > uJhoUow, the self- 
similar solution becomes "hollow:" There exists some fi- 
nite £i„ > such that the spatial region £ < £i„ is evacu- 
ated (p = 0). The self-similar solution describes the flow 
for £i„ < £ < 1 and is connected to the evacuated region 
£ < £m by a weak or tangential discontinuity, which lies 



at £ = £i„. In this case, U tends to 1 and C tends to 
as £ tends to £i„. 

Let us examine in more detail the uj < oJhoiiow solu- 
tions. It is straight forward to show that in the limit 
£ — > we have, to leading orders in £, 

[/(£)= 1/7 +c/ir, 
c(£)=Cor"/' + Cir/', 

G(£) = Gor-' + Gie2-2, 

P{0 = Po + PiC, (11) 



where 



7-1 



(12) 



Note that a > 2 for a; < 3/7 (and 7 > 1). Hence, the en- 
tropy (and temperature) diverge as £ tends to zero. This 
divergence is due to the fact that in the self-similar so- 
lution the shock becomes infinitely strong (infinite Mach 
number) as i? — > 0, which implies that the entropy of the 
gas shocked at i? tends to infinity. For any physical 
flow, where the energy is released within a finite radius, 
the entropy approaches a finite value as r ^ 0. This is 
the origin of the difficulties in choosing proper boundary 
conditions for the stability analysis (see § 13. 2|) . For the 
case 3/7 < < uJhnUov} the flow is conve ctively unstable, 
as indicated bv lRvu fc Vishniael l)199H) . 

3. PERTURBATION ANALYSIS 

We first derive the equations describing global pertur- 
bation evolution in § 13.11 and then discuss the proper 
boundary conditions in § 13.21 

3.1. Perturbation equations 

We use here the Eulerian perturbation approach, i.e. 
we define the perturbed quantities as the difference 
between the perturbed solution and the unperturbed 
one in the same spatial point. The derivation of the 
perturbation equation i s sim ilar to the one given in 
iRvu & VishniacI l|1987l 1199 ID . We therefore give here 
only the definitions and the main results. We define 



Su{r, e, LP, t) = u(r, e, ip, t) - Uo(r, t)f , 
Sp{r,0,ip,t)^p{r,0,ip,t) - po{r,t), 
Sp{r, 0, LP, t) =p{r, 6, (f, t) - po{r, t). 



(13) 



where u, p, and p are the velocity, pressure and 
density in the perturbed solution and uqY, pq, and 
po are the same quantities i n the unpertu rbed solu- 
tion. In this ana l ysis, a s in iR^fc VishniacI llT987i 
1991. ICherolien llT99fih IGoodmanI (1199(11) , and 
Sari. Waxman fc Shvarts I l|20n(lD . onlv "global" pertur- 
bations are considered, i.e, perturbatio ns that can be 
written in a separation-of-variables form l|Cox Ill980|) : 

Su{r, 0, ifi, t)=^R[SUr{0Yi,„ri9, ip)r - 5Ut[O^tYi^A0, v)]f[t), 
5p{r, 0, if, t)^BR'SGiOYi^^{0, ip)f{t), 

6p{r, 0, ip, t) = BR'RHPiOYi,miO, if) fit), (14) 
where 



Vt = 



80 



d_ 

sin dip 



(15) 



is the tangential component of the gradient. R{t) or sim- 
ply R is the unperturbed shock radius. The perturbed 



4 



shock radius, R(t, 6, ip), is given by 

R{t, 9, ip) - Rit) EE SR{t, 9, ^) = Rim^raiO, V)f{t), 

(16) 

and the dimensionless spatial coordinate, ^, is normahzed 
with respect to the perturbed shock radius, f = r/(i? + 
SR). 

Equations (|14|l and IjlGI) are the definitions of the quan- 
tities SUr, SUt, SP, 6G, and /. The quantity / mea- 
sure the ampHtude of the perturbation relative to the 
unperturbed values. If the function / is bounded for any 
time then the solution is stable, while if / is unbounded 
then the solution is unstable. The question of stability 
is therefore a question of solving for the function /. 

Linearizing the hydrodynamic equations around 
the unperturbed self - similar solution one obtains 
lIRvu fc Vishniac IIToSTI IChevalier lIToonl^ : 

SG{q-uj + 3U + £,U') + £,5G'{U - 1) + 5Ur{£.G' + 3G) 
+ Gi5Ul -l{l + 1)^G6Ut = 0, 
{5+q + 2U-l+ ^U')G^5Ur + {U~ 1)G^^SU^ 

{S +q + 2U- 1)G^SUt + {U- l)G(^6U^ + C^5P = 0, 
q — — "fq-pT- + (c/ — l)[—dF — "f-p;oG —oF 



P G ' "P ' G 
FC P' c 

+ l-^5G] + {^'l-^)mr^Q. 



p2 



where 



fR 



const. 



(17) 



(18) 



allows for variable separation. 

The value of q, defined in equation H18|l . is in gen- 
eral a complex number. The perturbation amplitudes 
are, in general, also complex, and the physical solution 
is obtained by taking the real part of the complex solu- 
tion. The time dependence of the perturbation ampli- 
tude, / oc i?^ oc t°"^, can be written as 

f oct^'"^'^ exp[ilm{s)\nt] with s = aq. (19) 

Thus, the real part of s describes a power law growth (or 
decay) of the perturbation amplitude, and the imaginary 
part describes oscillations of the perturbation amplitude 
with log(time). Positive (negative) values of Re(g) cor- 
respond to unstable (stable) perturbations. 

3.2. Boundary conditions 

The value of the parameter q is determined by the 
boundary conditions. At the shock front, the linearized 
Rankine-Hugoniot jump conditions are expressed as 

5G(l) = -^^w-G', 



SPil) 



7-1 
^ q-U', 



7 + 1 
2 

7 + 1 
2 

7 + 1 



For the case of 3 > w > tOhoUow, where the unperturbed 
solution is evacuated within < ^in, the proper (phys- 
ical) boundary condition at an interface between fluid 
and vacuu m is that the L agrangian pressure perturbation 
vanish (Goodman 1990). For the case oi u; < uihoiiow, 
■Rvu fc Vishniac (1987. 1991.1 required JP(0) = in or- 
der for the fluid not to undergo divergent perturbations 
in the Lagrangian sense at the origin. However, as ex- 
plained in § 12.21 the unperturbed solution shows non 
physical divergences at the origin, and the physical mean- 
ing of applying boundary conditions at ^ = is unclear. 

In order to overcome this problem, we use the follow- 
ing method. As explained in § 12.21 the entropy of any 
physical flow approaches a finite value as r ^ 0, while 
the entropy of the self-similar solution diverges at this 
limit. Thus, instead of directly analyzing the global sta- 
bility of the decelerating shock waves through an analy- 
sis of the Sedov-von Neumann- Taylor solutions, we an- 
alyze the stability of modified solutions which deviate 
from the Sedov-von Neumann- Taylor solutions in a re- 
gion near the origin, < ^ < S,o{t) with ^o(0 — > as 
i? — > cxj (where, e.g., the entropy of the modified solu- 
tion is finite). As we show below, the stability properties 
of the solutions are independent of the details of the so- 
lution at < ^ < ^o(^)i and are completely determined 
by the requirement that the perturbed solution remains 
self-similar at ^ > (t) ■ The stability of the Sedov-von 
Neumann- Taylor solutions is obtained by examining the 
solutions in the limit ^ 0. We show below that the de- 
rived growth rates converge to finite values in this limit. 

Since the solution we consider deviates from the self- 
similar solution at ^ < S,o{t), a discontinuity (of the en- 
tropy or of one of the derivatives of the velocity and 
pressure) occurs at (t) ■ Such a contact or weak discon- 
tinuity propagates along a characteristic of the flow, i.e. 
^o(^) must be a characteristic. Since G+ characteristics 
of the Sedov-von Neumann- Taylor solutions originating 
from any point £ < 1 ove rtake the shock front £ = 1 
(jWaxman Sz Shvarts Ill993|) . £o(t) must coincide with a 
Go characteristic of the Sedov-von Neumann- Taylor so- 
lution. The equation describing the time evolution of Go 
characteristics is 



(21) 



which may be written also as (jWaxman fc Shvarts 1199^ 



log go 
dlogR 



1. 



(22) 



For go Eq. (|22|l can be approximated using eq. (|ll|l 



dlogR 7 



l^£o^i?i-i/7_i"(i-i/7), (23) 



[2(q + 1) - c^] - P'. 



(20) 



Thus, at late times go tends to zero, i.e. the fractional 
size of the non-self-similar region tends to zero. 

The required additional boundary condition for the 
perturbed solution is obtained using the following ar- 
gument. The unperturbed solution is described by the 
self-similar solution in a region £,q{R) < g = r/R < 1. 
rc{t) = goi? is the radius of a sphere across which there 
is no mass flux, and within which the solution deviates 
from the self-similar solution. In the perturbed solution, 
the surface across which there is no mass flux is deformed 
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and its radial distance from the origin depends on direc- 
tion, fciO,(f>,t). The evolution of this surface with time 
is given by 

(IT ' " ~ • ~ ~ 

" 0, 0, t) = B^oUi^o) + fmoSUr{^o)Yim{0, ^) 



dt 



+ io5R-^{RioU{io)), 



(24) 



where = fcl{R + 5R), Ur is the perturbed radial ve- 
locity (the transverse flow leads only to second order cor- 
rections to the evolution of fc), and the last two terms 
on the right hand side are the Lagrangian variation in 
ua{£,oR,t) — R^oU{^a). As for the unperturbed solu- 
tion, we assume that a perturbed solution may be con- 
structed, such that it deviates from self-similarity only 
within r{9, (j), t) < fc{9, (p, t). In order to obtain a global 
self-similar perturbed solution, defined over the same ^ 
range as the unperturbed solution, we require ^o(^) to 
coincide with £,o{R), i.e. [see eq. (|22|) ] 



log Co 
dlogii 



or [see eq. (|^ ] 



d{R 



dt 



(25) 



(26) 



Comparing equation H24|l with equation H2t)|) . and using 
^0 = £,0j we obtain 

SUriCo)=qU{^o)-^oU'i^o)- (27) 

Eq. (|27ll provides the additional constraint required for 
determining q: Integration of equations H17() starting at 
the shock front using the shock boundary conditions H2()|l 
will lead to solutions satisfying the condition H27|l only 
for some particular values of q. 

The procedure described above determines 
g(/, 7, w, Co(0)- The time dependence of implies 
a time dependence of q, while the derivation of eqs. lfT7|l 
is valid under the assumption that q is time independent. 
The corrections to eqs. (I17II due to the time dependence 
of q may be neglected if perturbations evolve much faster 
than the rate at which q changes, i.e. for \t/s\ ^ |s/s| 
which may be written as 



|(7|»(l-;7(eo))| 



dlog(g) 



(28) 



'rflog(eo)' 

As we show in the following section, s = qa converges to 
a finite value in the limit 0. Thus, the value of s 

in this limit provides the growth rate of global perturba- 
tions to the self-similar solutions. Moreover, the value of 
s obtained for finite gives the growth rate for global 
perturbations to solutions which deviate from the self- 
similar solution at ^ < provided the condition given 
in eq. H28|l is satisfied. 

4. NUMERICAL RESULTS 

We present in this section the results of numerical cal- 
culations of s{-f, uj, I, £_o). Eqs. ((T7|l are integrated start- 
ing at the shock front, with boundary conditions (|20|) . 
using standard error control Runge-Kutta integration 
scheme. The value of s is obtained, using Newton- 
Raphson iterations, by requiring the boundary condi- 
tion H27|) to hold. In addition to presenting results as 




Fig. 1. — Perturbation growth rate, s, as function of wave number 
for various values of mo , the fraction of self-similar solu tion mass 
contained in the self-similar part of the flow [see eq. II29I 1 , for uj = 
0,7 = 1.5 . The lines deno t ed "R fcV" show the results of the 
analysis of IHvu fc Vishniac I 119871) . The perturbation amplitude 
evolves as / oc t^'^^"^ exp(ilm(s) Int). 




Fig. 2. — Perturbation growth rate, s, as function of wave number 
for various values of mo , the fraction of self-similar solution mass 
contained in the self-similar part of the flow, for u! = 0,-f = 1.1. 
The lines denoted " Rfc V" show the results of the analysis of 
Evu & Vish niaci 1198711 . The perturbation amplitude evolves as 
/ oc i^^^f") exp(ilm(s) Int). 



function of ^q, we also present results as a function of 
'^o('?o)j where ttiq is the fraction of self-similar solution 
mass enclosed within the self-similar region, 



mo 



(29) 







Some values of for different rng's and 7's in the uj 
case are shown in tabled 

The results for ui — 0, i.e. the frequency spectra of 
the real and imaginary parts of s, are shown in figures ^ 
and |2] for j — 1.5 and 7 = 1.1 respectively. The spa- 
tial dependence of the perturbations to the self-similar 
profiles of flow variables are shown for these cases in fig- 
ures 13 to |S1 The convergence of s to finite values in the 
limit jTiQ — > 1 justifies the use of the method described 
in the previous section (see last paragraph of ii l3.2|l . The 
self-similar, asymptotic solutions are stable for large 7, 
and unstable for small 7. The various stability regions 
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TABLE 1 

Some values of go for different mo's and 7's in the o; = case 



mo 


7 = 5/3 


7 = 1.5 


7 = 1.4 


7 = 1.3 


7 = 1.2 


7 = 1.1 


7 = 1.06 


0.999 


0.4312 


0.4965 


0.5494 


0.6170 


0.7061 


0.8273 


0.8886 


0.99 


0.5850 


0.6411 


0.6839 


0.7366 


0.8026 


0.8874 


0.9280 


0.95 


0.7207 


0.7629 


0.7940 


0.8309 


0.8756 


0.9306 


0.9564 


0.9 


0.7855 


0.8196 


0.8443 


0.8733 


0.9076 


0.9490 


0.9681 


0.85 


0.8245 


0.8533 


0.8740 


0.8980 


0.9261 


0.9595 


0.9747 


0.8 


0.8523 


0.8772 


0.8949 


0.9152 


0.9388 


0.9667 


0.9792 


0.75 


0.8738 


0.8955 


0.9108 


0.9283 


0.9485 


0.9721 


0.9826 


0.7 


0.8913 


0.9103 


0.9237 


0.9388 


0.9562 


0.9763 


0.9853 


0.65 


0.9059 


0.9227 


0.9343 


0.9475 


0.9625 


0.9798 


0.9875 


0.6 


0.9185 


0.9332 


0.9434 


0.9548 


0.9678 


0.9827 


0.9893 



(0=0,->tl.5 , /: 



(0=0, 1^1.1 , 




(0=0, -pLs , ;=io 



(D=0,ltl.l ,1=40 



Fig. 3. — The real and the imaginary parts of 5G as function of 
1 — whe re m (g) is the fraction of self-similar solution mass in 

[0,g] [see eq. 1^ 1. for the = 0,7 = 1.5, Z = 10 case {(a)-{b)) and 
for the oj = 0,7 = 1.1, Z = 40 case ((c)-(d)). The physical solution 
oscillates between the real and the imaginary parts. 
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Fig. 5. — The real and the imaginary parts of SUr as function of 
1 — m(g), whe re m (g) is the fraction of self-similar solution mass in 
[0,g] [see eq. 1^ 1. for the = 0,7 = 1.5,1 = 10 case {(a)-{b)) and 
for the w = 0,7 = 1.1, Z = 40 case ({c)-(d)). The physical solution 
oscillates between the real and the imaginary parts. 
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Fig. 4. — The real and the imaginary parts of 5P as function of 
1 — m(g), whe re is the fraction of self-similar solution mass in 

[0,g] [see eq. 1^ 1. for the = 0,7 = 1.5,/ = 10 case {(a)-{b)) and 
for the w = 0, 7 = 1.1, Z = 40 case ((c)-(d)). The physical solution 
oscillates between the real and the imaginary parts. 
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Fig. 6. — The real and the imaginary parts of SUt as function of 
1 ~ m(^), whe re m (g) is the fraction of self-similar solution mass in 
[0,5] [see eq. 1^ 1. for the (.j = 0,7 = 1.5,1 = 10 case {(a)-{b)) and 
for the a; = 0,7 = 1.1,/ = 40 case ({c)-(d)). The physical solution 
oscillates between the real and the imaginary parts. 



in the [^,7] plane are shown in fig. [T] For all values of 
7, the perturbation amplitude shows oscillatory behavior 
(Im(s) ^ 0) at all but the lowest {I = 1) wave numbers. 



These oscillations are dumped (Re(s) < 0) in the stable 
solutions, and amplified (Re(s) > 0) over a wide range of 
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Fig. 7. — Stability regions in the [oj, 7] plane: Solid thick line- 
the asymptotic stability line (mo — > 1); Dotted line— uj = 8/7, 
above which the solutions are convectively unstable. The dashed 
thick line, uj = i^hollomil) shows the border of the region (top right 
part of the plane) where the self-similar solutions are hollow, and 
where our analysis does not apply. 



wave numbers (10 < I < 100) in the unstable solutions. 

The figures demonstrate that a value of /, / = Iq, 
exists (Zo = 2 for [w = 0,7 = 1.5] and lo = 10 for 
[uj — 0,7 = 1.1]) such that for I > Iq the solutions 
converge to the Ryu & Vishniac solutions as mo tends 
to 1. The reason for this convergence is explained in 
appendix § In brief, there is some component of 
the perturbed solution that diverges at the origin. In 
this situation, any finite boundary condition forces this 
component to vanish, thus yielding the correct solution. 
For I < lo, the results of Ryu & Vishniac differ from 
ours. I t is interesting to note here that it was proven by 
IGurevic h & Rum yanste v (1970), using arguments valid 
for Z — > and any [lo,j], that two stable modes of per- 
turbation exist in this limit, with s = and s = — 1 (for 
details, see appendix §|b1i. Our method of solution re- 
produces this result, in contrast with the method of Ryu 
& Vishniac, which does not provide the correct value of 
q for small wave numbers. 

As explained at the end of S I3.2I the value of s obtained 
for mo < 1 gives the growth rate for global perturbations 
to solutions which deviate from the self-similar solution 
at ^ < ^0(^0)7 provided the condition given in eq. 128|) is 
satisfied. Figure |S1 shows a magnification of figure|7|with 
contours (dashed lines) indicating values of mo for which 
equality is obtained between the left hand side and the 
right hand side of eq. H28(l . for the wave number 
at which a maximum in |s| is obtained. The growth 
rates s obtained by our analysis are accurate for mo val- 
ues larger than those at which equality is obtained. For 
small 7 values, e.g. 7 < 1.06 for cj = 0, our analysis pro- 
vides accurate s values for arbitrarily small mo. Figure |H1 
demonstrates that, based on the current analysis, one 
can not infer the existence of instability for mo < 1 at 
[a;,7] regions where the asymptotic, mo 1, solutions 
are stable. However, we find that the flow is in gen- 
eral less stable for smaller mo. In the region where the 
asymptotic solutions, mo 1, are stable, we find slower 
decay rates of the perturbation amplitudes for mo < 1, 
compared to the decay rates obtained for the asymptotic 
solutions. In the region where the asymptotic solutions 
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III I 
III I 

'1 1.1 1.2 1.3 1.4 1.5 
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Fig. 8. — Regions of validity: The growth rates derived by our 
analysis are valid for solutions that coincide with the self-similar 
solution at ^0 < € < Ij provided that the fraction of self-si mila r 
solution mass contained in the [§o, 1] region, mo(5o) [see eq. 1291 1. 
is larger than the value indicated by the dashed lines. The thin 
solid line shows the boundary of the region where the gr owth rates 
derived by our analysis are valid for any value of mo [eq. I28i holds 
for any mo] . The thick solid line shows the boundary of the stability 
region of the asymptotic solutions, mo — » 1, and dotted lines show 
the marginal stability lines for different mo values, obtained by 
our analysis. A comparison of these lines with the dashed lines 
shows that instability at mo < 1 is obtained in the region of stable 
asymptotic solutions only for mg values lower than those for which 
our analysis provides accurate values of the growth rates. This 
implies that, based on the current analysis, one can not infer the 
existence of instability for mo < 1 at regions where the asymptotic, 
mo — > 1, solutions are stable. 



are unstable, we find growth rates which are larger than 
those obtained for the unstable asymptotic solutions. 

5. DISCUSSION 

We have examined the stability of decelerating shocks 
expanding into ideal gas with density profile decreasing 
with distance from the origin according to p = Kr^^ , 
for the case where the solutions are not hollow, i.e., lo < 
(7 — 7)7(7 + 1). We have shown that direct examination 
of the stability of Sedov-von Neumann- Taylor solutions is 
not possible due to the divergence of the entropy near the 
origin (see §[221and §EIl. Instead, we have analyzed the 
stability of modified solutions, that coincide with the self- 
similar solutions at < ^0 < C = '''/Ri't) < 1; ^^^d deviate 
from the self-similar solutions at ^ < ^o- We have shown 
that the global stability of the flow is independent of the 
details of the solution in the region where it deviates from 
self-similarity, and obtained the growth rates s{l,^,uj) 
[see eq. ()19(l ] of global modes of perturbations by taking 
the limit ^0 ~* 0- Several examples of the dependence of 
s(Z, 7, bj) on I and of the spatial structure of global modes 
of perturbations are given in figures ^ to |B1 The regions 
of stability in the [0^,7] plane are outlined in figured 

Using our new method of analysis, we have demon- 
strated that while the growth rates of global modes de- 
rived by previous analyses are accurate in the large wave 
number (small wavelength) limit, they differ significantly 
from the correct values at low wave numbers (see, e.g. 
figures n and |2J). The reasons for the convergence of 
previous analyses (or any other analysis which require 
finite inner boundary condition) to the correct values of 
s(Z,7,w) at large /, and for their failure at small /, are 
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-2 -1 1 2 
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Fig. 9. — The real part of the perturbation growth rate, s, as 
function of wave n umber for a; = 0, 7 = 1.06, compared with the 
values measured bv lGrun et al. 1 119911) . The measured values are 
higher than the growth rates obtained for the asymptotic, mo — > 1, 
self-similar solutions. However, as explained in the discussion, the 
flow in the experiment reaches only mo < 0.87. 



explained in appendix A. 

The growth rates obtained by our method for finite 
s(Z,7,w;^o)j describe the growth rates of global pertur- 
bations for solutions that deviate from the self-similar 
solutions at ^ < £,q, provided the constraint of Eq. H28|l 
is satisfied. Figure |H1 shows the values of mo{uJ,j), the 
fraction of self-similar solution mass contained in the self- 
similar part of the flow, ^0 < C < 1 [see eq. (I29ll ]. above 
which this condition is satisfied. For small 7 values, e.g. 
7 < 1.06 for w = 0, our analysis provides accurate s 
values for arbitrarily small mg. 

Based on the current analysis, one can not infer the ex- 
istence of instability of the pre-asymptotic flow, i.e. for 
TTiQ < 1, at [0^,7] regions where the asymptotic, mo — > 1, 
solutions are stable (see figurelH)). However, we find that 
the pre-asymptotic flow is in general less stable than the 
asymptotic flow. In the stable region of the [uj,j] plane 
(figure IHl), we find (see figure Q]) slower decay rates of the 
perturbation amplitudes for toq < 1 (Co > 0), compared 
to the decay rates obtained for the asymptotic solutions, 
mo — *■ 1 (Co 0). In the unstable region, we find growth 
rates which are larger than those obtained for the unsta- 
ble asymptotic solutions (see figure This implies that 
as the shock wave evolves and approaches self-similarity, 
the perturbation growth rates are also evolving, and are 
determined at any given time by the appropriate value 
of mo. For unstable shock waves, the over all growth of 



perturbations depends on the evolution of mo at early 
times. 

I ■ the experiment of iGrun et al. I l|199H) . the ambi- 
ent gas adiabatic index was 1.06 ± 0.02, and the max- 
imum growth rate was s = 1.6 at kR = 10 (we will 
assume kR ~ I though the wave numbers in the exper- 
iment were obtained from the projection onto a plane 
of the edge of an unstable sphere, thus kR is not iden- 
tical to I). The numerical and the experimental results 
for w = and 7 = 1.06 are shown in figure © (Note, 
that our method is valid in this case for any mo, see fig- 
ure The maximal growth rate of perturbations to the 
asymptotic solution is s = 1 at Z = 60. The growth rate 
is larger, however, for mo < 0.7. W e can estimat e mn i n 
the experiment from figure 3(a) of IGrun et al. I (|1991"). 
The radius beyond which the shock expands according 
to the self-similar scaling is RmiUai — 7mm, and the 
perturbations keep growing until t = 200ns, where the 
shock radius is R final — ^'^mm. Thus, toq is given by 
m^{R) = 1 — {Rinitiai / R)^ , and the maximum value of 
mo is mg ~ 0.87. This implies that the experiment does 
not measure the asymptotic (toq — > 1) value of s, but 
rather s(too) with mo < 0.87. This may account for the 
large growth rates determined by the experiment. 

The most unstable modes predicted by our linear anal- 
ysis are still characterized by larger wave numbers, I > 
50, than obtained by Grun et al. (1991). However, since 
the growth rates predicted by our analysis for large wave 
numbers are large, we expect the evolution of perturba- 
tions to become non-linear at an early stage. The non- 
linear interaction at large wave numbers may lead to per- 
turbation amplitudes at smaller wave numbers, I ~ 10, 
which are larger than predicted by linear analysis. Our 
predictions may be tested in future experiments, by con- 
trolling mo- For example, lower gas densities or higher 
ablated material densities will make mo smaller, leading 
to higher perturbation growth rates. 

Finally, we note that for stable asymptotic shocks, e.g. 
the = 0,7 = 5/3] case relevant for the " Sedov-phase" 
of supernovae shocks, our analysis results do not differ 
significantly from from those of Ryu & Vishniac. How- 
ever, since in general we find that the flow is less stable at 
small mo, it is possible that instabilities do arise during 
the pre-asymptotic stage of shock evolution. This should 
be checked by numerical simulations. 

DK thanks Re'em Sari for discussions that triggered 
his interest in this problem. This research was partially 
supported by AEG and MINERVA grants. 



APPENDIX 

A. THE ASYMPTOTIC BEHAVIOR FOR LARGE Mo VALUES 

As indicated in §^ our solution for the growth rate s tends to the Ryu & Vishniac solution as mo tends to 1, for 
sufficiently large Z, / > Iq. We explain here the reason for this convergence, and the reason for the deviation of our 
results from those of Ryu & Vishniac for I < Iq- Moreover, we will prove that any finite inner boundary condition will 
yield the correct solution for sufficiently large I. 

The perturbation eqs. H17|l can be simplified in the limit ^ — * by using f following iRvu Sz Vishnia'cl l)1987j) ) the 
asymptotic behavior, Eq. Hll|) . of the self-similar solution in this limit. Assimring that in this limit the perturbations 
behave as 

^^ = ^^0^", SUr = SUro^''-% SUt^SUto^''-'', 60 = 500^"-^ (Al) 
a set of four linear equations is obtained for 6Pq, SUro, SUtoi a-nd 6Gq. In order for a nontrivial solution to exist, the 
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determinant of the coefficient matrix of these equations must vanish, 



All 


A12 


Ai3 





A21 


A22 





A2i 








A33 


Asi 


Ail 


Ai2 









where 



An: 

A12: 
A13: 
A21: 
A22 

A2i- 

A33: 

A34: 
A41: 

Ai2-- 



0, 



+ 3/7+(d-2)(l/7-l), 



(A2) 



1 



1 



(l/7-l)(d-a), 
(l/7-l)(ci-a), 



-q — uj 

--{d+1), 
---1(1 + 1), 
---aPi/Go, 
--6 + q + 2/-/ 
-d, 

--6 + q + 2/-f 
-A, 

---jq+{j-l){d-a), 

= -7(a-2). (A3) 

Eg. IA2l is therefore a fourth order polynomial equation for d. If four different solutions d exist, none of the powers in 
Eq. lAll is zero, then any solution of the perturbation equations is given, near the origin, by a linear combination of 
the power-law solutions corresponding to different values of d. For I > Iq, f our d ifferent solutions to the polynomial 
eq. for d exist, with only one negative d — a value (the smallest power in Eg. IA1|) . Requiring the pressure to be finite 
at the origin is equivalent in this case to requiring that the solution be a combination of the 3 approximate solutions 
corresponding to positive d — a values. This uniquely determines the value of q. Thus, the method used by Ryu 
& Vishniac, seeking solutions that do not diverge near the origin by requiring finite boundary condition, yields the 
correct value of q for I > Iq. 

In the region I < Iq, more than one solution with negative d — a exists. In this case, requiring the solution not to 
diverge at the origin does not yield the correct solution. It is for this reason that the numerical method used by Ryu 
& Vishniac does not converge to the correct solution of q for I < Iq. I n con trast, our numerical method is free of such 
difficulties. Note, that in the limit / 0, one o f the solutions to Eg. IA2I is d = 0, which implies the existence of an 
additional solution, not of the form given by Eg. I All This additional solution is given analytically in appendix S IbI 

B. THE BEHAVIOR FOR SMALL WAVE NUMBERS 

In the self-similar solution, any flow variable w{r,t) (where w may stand, e.g., for p,u,c) is given in the form 
w{r,t) = R^^W{^) where ^ = r/R. For the problem under consideration, R is given in the form R — A{t — to)". 
The parameters ei and 62 and the function W depend only on the parameters 7 and u (which determine also the value 
of a). Thus, for given j,uj the solution depends on only two parameters, A and <o- 

An approximation to the perturbed solution in the limit of small wave numbers, / —>■ 0, may be obtained using the 
following argument ijGurevich &: Rumvanste'v1ll970fl . For small wave numbers, large wave lengths, we may assume 
that the flow along any direction 0, (j) evolves with time as if it were part of a spherically symmetric solution. This 
approximation is valid as long as the time scale for the evolution of the perturbation is shorter than the time it takes 
sound waves to propagate tangentially over a distance ^ R/l over which the deviation from isotropy is appreciable. 
This approximation gives, of course, the exact solution for I = 0, and may be expected to provide an approximate 
description of the behavior at I —>■ 0. Since the spherically symmetric solution depends on only two parameters, A and 
to, the perturbed solution is given in this approximation in the form 

wir,e,ip)^R;'R;-W{r/Rp), (Bl) 

where 

Rp^R{t;Ai0,if),toie,ip)). (B2) 

Here, it is assumed that A{9,ip) and to(0,ip) differ from their values in the spherical solution, Ag and ips, only 
infinitesimally, and that they are slowly varying functions of {0, ip). 
Using spherical harmonics decomposition, 

Ai0,ip)=A + SAYimie,ip), 

to{0,v)^to+5toYi„{e,(p), (B3) 

one obtains 

6R{t, 9, ip) = [—SA + —Sto]Yi„,{0, ifi), 



5w{r,0,ip)^R''^R''^[ei-^W' + 62-. ]SR 

RSR 



(B4) 
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for the deviatfon from the spherical solutfon. 

From equation (|B4(I it is clear that there are two perturbation modes: one is associated with changes in to, and the 
other with changes in A. Since for R{t) = A{t — to)" we have SR{t) — R{SA/A + Sto/{t — to)), s = is obtained for 
perturbations associated with changes in A, and s = — 1 is obtained for perturbations associated with changes in to- 

The radial dependence of the perturbations is given by 



for to associated perturbations, and by 



SGiO = [-uj-CG'{0/GmG{0, 
6UriO = [S-l-CU'{0/umu{0, 

sp{o=[-cj +26 - iP'io/m]m 



(B5) 



(B6) 



for A associated perturbations. 

The tangential velocity perturbation may be obtained by integrating the third equation of H17|l . Since it is propor- 
tional to VTYim{0,y^), the tangential velocity perturbation tends to zero as I tends to zero. It is straight forwar d to 
show that for every lu and 7 equation (jB5|) is a solution of equations H17|) with I = and s — —1, and that equation (jB6p 
is a solution of equations H17() with I ~ and s = 0. This can be most easily inferred from equation lO, equation (jS)), 
and energy conversation, 

^h~i)uHi-u) 
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